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We present results from simulations of two flavors of dynamical overlap fermions on 8* lattices 
at three values of the sea quark mass and a lattice spacing of about 0.16 fm. We measure the 
topological susceptibility and the chiral condensate. A comparison of the low-lying spectrum of the 
overlap operator with predictions from random matrix theory is made. To demonstrate the effect of 
the dynamical fermions, we compare meson two-point functions with quenched results. Algorithmic 
improvements over a previous publication and the performance of the algorithm are discussed. 
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I. INTRODUCTION 

On i . . 

, Chiral symmetry is a fundamental part of the theory of strong interactions and therefore should be respected when 
' putting QCD on the lattice. It was realized more than two decades ago, that this can be done with Dirac operators 
. D which (at zero quark mass) obey the Ginsparg- Wilson equation 1] 

c : a 
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where a is the lattice spacing and Rq the radius of the Ginsparg- Wilson circle. A decade ago, this constraint was 
' realized by overlap fermions f2^)'3|, fixed-point fermions 4] and domain wall fermions 5, 6]. Among those, only domain 
I wall fermions have been used for some time in simulations on the lattice which include the fermionic determinant. 
^ Only recently the first steps toward dynamical simulations using overlap fermions have been taken 0, 0, IS liol ■ 
T-H , Because of the high cost of applying the Dirac operator these are still limited to a small volume. In this paper we 
' present results from simulations with two flavors of dynamical overlap fermions in a small box and at a large lattice 
' spacing. We measure the topological susceptibility and, by comparing the low-lying eigenvalue distribution to random 
matrix theory, the chiral condensate. We also fix the lattice spacing and look for dynamical effects in meson two-point 
functions. The measurement of the susceptibility is greatly facilitated as compared to, e.g., simulations using highly 
improved staggered fermions 'l2]. With the overlap operator, we can use the same Dirac operator in the simulation 
"j*-^ ' and in the definition of the topological charge via the index theorem. The topological charge defined in this way thus 

, directly influences the weight of a conflguration. 
' I ■ A more technical consequence of that is the discontinuity of the fermionic determinant as a function of the gauge 
CIh: variables. The surfaces in gauge fleld space at which the fermionic action is discontinuous coincide with the change in 
^ ' topology defined by the index of the Dirac operator. Although this does not pose a problem in principle, in practice 
""T^ ■ changing the topological sector turns out to be a significant problem. As recent simulations with highly improved 
I> . stag gered fermions show, this problem itself, however, is not restricted to chiral formulations of QCD on the lattice 

; S _ 

, In a recent paper [l^l, we described the setup of our simulations of two degenerate flavors of overlap fermions. We 
' studied the impact of fat (stout !l^) links and multiple pseudo-fermion fields. Some improvements to the algorithm 
will be reported in the following. The simulations presented in this paper used very limited computer resources, i.e. 
half a year on an array of 12 3.2 Ghz Pentium-IVE's. With that we present some results on 8* lattices and a relatively 
coarse lattice spacing of around a = 0.16 fm and a quark mass down to about 35 MeV. 

This paper is organized as follows: We first review the algorithm, describe the improvements over our previous 
publication and give the parameters of our simulation. Then in Sec. IIVI we attempt to set the lattice spacing by 
determining the Sommer parameter tq and proceed with the extraction of the topological susceptibility in Section 
In Section IVII we look at meson two-point functions. By comparing to quenched results on matched lattices, we 
demonstrate the impact of dynamical fermions on the scalar correlator. Finally, we compare the low-lying spectrum 
of the overlap operator with the predictions from random matrix theory. A note in caution: for all these investigations, 
the volume which we are simulating is too small. Future simulations will improve on this. Here we want to demonstrate 
that simulations with dynamical overlap fermions are possible and give results which match our expectations of full 
QCD in a small box. 
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II. DEFINITIONS AND ALGORITHM 



Let us fix the conventions and describe the algorithm together with the improvements since our_previous pubhca- 
tion [lol |. to which we refer the reader for more details. The massive overlap operator is given by 2,'S] 

7TL 

Dov{m) = (i?o - y) [1 + l^<K-Ro))] + m (2) 

with e{h) = h/\/h? the sign function of the Hermitian kernel operator h = 75^ which is taken at the negative mass 
Rq. Rq is the radius of the Ginsparg- Wilson circle. We are using a planar kernel Dirac operator d with nearest and 
next-to-nearest ("\/2") interactions. For details see Ref. IQ]. The sign function is computed using the Zolotarev 
approximation with an exact treatment of the low-lying eigenmodes |A) of h{—Ro) 

eihi-Ro)) = ^sign(A.)|A.)(A.| « Usign(A.)|A)(A| -^(1 - |A.)(A.|) . (3) 

2 l—l J l — l 

For our simulation we use the Hybrid Monte Carlo (HMC) algorithm as modified for overlap fermions by Ref. 
0. The effective action is given by 



=1 



with H (m) — Dovim)^ Dovim) the square of the hermitian overlap Dirac operator. Sg[U] is the gluonic action. The 
are the Upf pseudo-fermion fields used to include the contribution of the fermion determinant. The use of several 
of these fields has been suggested in It improves the stochastic estimate of the determinant. We studied its 

effects extensively in Ref. JOJ . 

With the HMC algorithm, an ensemble distributed according to this effective action is generated by updating the 
gauge fields using molecular dynamics trajectories with a final accept-reject step. At the beginning of the trajectory 
one chooses momenta tt conjugate to the gauge fields and also refreshes the pseudo-fermions. One then integrates 
the resulting equations of motion (treating SeS as the potential) numerically in some fictitious simulation time r. 
They result from the requirement that the total 'energy' TC = ir'^ /2 + Scs is conserved. At the end one applies an 
accept/reject step with acceptance probability P — min[l, exp(— A7i)] which corrects for the errors in the numerical 
integration. 

An important difference between conventional fermions and overlap fermions is that the effective action S'cff is 
discontinuous for the latter. The discontinuity has its origin in the sign function in the definition of the Dirac 
operator Eqs. Q and (|3J|. The fermionic action is discontinuous at the surfaces where the kernel operator h{—Ro) 
has a zero eigenvalue. According to the index theorem these are also the surfaces where the topological charge 
(as seen by the fermions) changes. Ref. Q gives the prescription for how to deal with this situation. The molecular 
dynamics evolution can be thought of as resembling that of a classical particle in the presence of a potential step. If 
the step in the action is too big for the particle to get across it, the particle is reflected, i.e. the momentum component 
normal to the zero eigenvalue surface is reversed. On the other hand, if the normal component is large enough to 
change topological sector, the normal component is reduced such that energy is conserved. Following Ref. ^ this is 
called a refraction. The momentum is then altered according to 



_ I -N {N\Tr) + N sign(7V|7r) y/{N\Tr)^ - 2ASf ii (N\tt)'^ > 2ASf 

^ j-2iV(iV|7r) if (iV|7r)2 < 2A5/ ^ ' 

with N the vector normal to the zero eigenvalue surface, tt the momentum and ASf the discontinuity in the fermionic 
action. To monitor whether an eigenvalue has changed sign, one thus has to compute some number rioig of the 
lowest eigenmodes of h(—Ro) and see whether the eigenvalue of any of them changes sign. The matching is done by 
computing the scalar products of the low modes before and after a step. Since the eigenmodes are needed anyway to 
precondition the construction of the sign function, there is virtually no overhead associated with this test. 

Note that this part of the algorithm potentially scales with the square of the volume: The cost of determining the 
height of the step is at least proportional to the volume. The number of times this procedure has to be executed can 
be assumed to be proportional to the density of eigenmodes of the kernel operator at the origin, which in turn might 
be proportional to the volume. It is therefore pivotal to keep the cost of this step as low as possible. 
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A major improvement in the algorithm is the way in which we compute the height of the step. In our previous 
pubhcation, we ran two conjugate gradients to compute 



AS' 



N-l 

E 

1=1 



(6) 



where the difference is taken between the fermionic actions for which only the sign of the lowest eigenmode (the one 
which becomes zero on the surface) is changed without changing any of the other modes or the gauge configuration. 
This was very expensive since we have to decide frequently whether to refract or reflect. Because the square of the 
hermitian Dirac operator in one chiral sector cr = ±1 is of the form 



Hlim)^2{Rl-—)P, 



l + f7^e(A0|A.)(A, 



(7) 



with |A) the eigenvectors of h{—Ro) and Pa = ^(1 + crjs) the projector on the chiral sector, the sign change of the 
crossing mode amounts to the change 



)Pa\Xo){\o\Pa 



(8) 



with |Ao) the zero mode and the sign being minus the product of the sign associated with the chiral sector and the 
sign of the eigenmode before the step. Thus, we can use the Sherman-Morrison formula to compute the height of 
the step. The result is 



A 



1 



H„{mf 



{ARl - m2) 



1 ± {ml - m-^){MP„Ha\'m)Pa\\o) 
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Ha{m)- 



rP.|Ao)|' 
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This has the additional advantage that instead of running two conjugate gradients for each pseudo-fermion field, we 
only have to invert once, using the eigenmode which changes sign as a source. Because the height of the step is directly 
computed, one also has better control over the accuracy — compared to taking the difference of two approximate 
quantities. 

One can further exploit Eq. by realizing that the cost the of the inversion of (m) depends on the chirality of 
the source and the topological sector in which one is inverting. One can compute the height of the step from either 
side of the surface. In the chiral sector with zero-modes, it is cheaper to invert in the topological sector of lower 
charge. However, since the zero- modes push the spectrum of the other modes up, the conditioning number of H^{m) 
in the sector without zero-modes is lower on the side of the surface with higher topology. A second advantage of the 
use of this formula is that we can monitor the step-height during the CG iterations and terminate the iteration when 
it has become clear that we are going to reflect. 

Finally, let us report a small improvement in the computation of the fermion force SSf[U]. The derivative of the 
Zolotarev part of the approximation to the sign function Eq. has been given in many places. For each pseudo- 
fermion field and each order in the rational approximation the formula has a term 1/ (ft,^ — c), the derivative of which 



- — {hSh + 5hh)—^ 



(10) 



One thus has to invert the kernel action, which can be done simultaneously for all shifts using a multi-mass algorithm. 
The computation of 5h follows via standard methods from Ref. j23|. However, due to the many shifts in the Zolotarev 
approximation and several pseudo-fermion fields, this part is a non-negligible contribution to the total cost of the 
simulation. 

More difficult is the derivative of the projector term P\ = |A)(A|. This derivative is basically given by first-order 
perturbation theory (See Ref. [21^) 



5Px 



^ {l-P^)5hP), + P),Sh+-^{l-Px). 



\-h 



(11) 



Because (A — /i) is singular, its inversion is problematic even though the contribution of the eigenmode with eigenvalue 
A is projected out of the source: since both A and |A) are only approximately known, one faces a "zero divided 
by zero" problem. In our previous publication, we therefore shifted the pole position, performed the inversion, and 
interpolated our result. This turned out to be insufficiently stable. Now we use a Chebychev approximation to the 
inverse of h'^ — A^ in the range such that the inverse is precise outside the known eigenvalues of /i(— i?o)- The advantage 
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TABLE I: Overview over our collected data at /3 = 7.2 with several streams per quark mass. We give the number of pseudo- 
fermion flelds Upf used, the number of trajectories, the acceptance rate and the total number of reflections and refractions. In 
our analysis we typicaUy discard the first 100 configurations per stream. The statistics includes all (accepted and rejected) 
trajectories. Note that the number of refractions is higher than the number of effective changes in topologies because of 
tunneling back and rejections. 



is that this approximation is finite at ^ = 0. The problem is thus reduced to a contribution from the eigenvalue 
A mode which is zero times some finite number given by the required accuracy and the range in which one computes 
the eigenvalues explicitly. 

Eq. [71 provides us with a tantalizing result, which we do not know how to apply in the context of HMC: an exact 
formula for the ratio of the fermion determinants on either side of the topology-changing boundary: 

ii|||i = l±,«S-„.^)(A.|P,^P,|A.). ,12) 

HMC does not ever use the exact determinant as part of the simulation. Instead, it generates configurations whose 
statistical weight is controlled by the effective action Eq. ^ In an algorithm which does approximate the determinant 
directly, like the Rq algorithm, it seems straightforward; one would just use the logarithm of Eq. E| as the step. 
However, we are unwilling to abandon HMC for two reasons: First, we feel that we benefit substantially from the fact 
that HMC is exact and therefore do not want to run an R type algorithm which cannot be made exact in an obvious 
way. Second, we gain considerable speed using HMC over an R algorithm because in HMC we can use a previous 
solution to an inverse of H{m)^ to begin the computation of the new force. This point deserves further investigation. 



III. SIMULATION PARAMETERS AND PERFORMANCE OF THE ALGORITHM 

We simulate on 8'' lattices at one value of the gauge coupling (3 — 7.2 which we chose to be roughly at the Nt = 6 
phase transition (for and overview over the lattice spacings see Table iHlj) . We use the Liischer-Weisz gauge action 12^ 
with the tadpole improved coefficients of Ref . 23] . Instead of determining the fourth root of the plaquette expectation 
value uq — {{Upi) self-consistently, we set it to 0.86 for all our runs as we did in our previous pubhcation. 

Our kernel operator d is constructed from gauge links to which two levels of isotropic stout blocking fl5| have been 
applied. The blocking parameter p is set to 0.15. 

We report on simulations at three values of the bare sea quark mass am^ = 0.03, 0.05 and 0.1. Based on measured 
lattice spacings from the Sommer parameter and the perturbative calculation of matching factors reported in the 
Appendix, we believe that these values correspond to MS quark masses of about 35, 55 and 100 MeV. An overview of 
our collected statistics is given in Table The trajectories all have length one, divided in 20 elementary time steps. 
We use a Sexton- Weingarten f24| integration scheme in order apply a smaller time step for the gauge field integration. 
We perform 12 applications of the gauge force and gauge field update per elementary time step. Again, see Ref. ^3 
for details. 

In order to monitor whether an eigenvalue has changed sign, we compute the lowest 8 eigenmodes of ft.(— i?o) in 
each step. These are also used to precondition the construction of the sign function. 
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FIG. 1: The stochastic estimate of the height of the step compared to the actual change in the logarithm of the determinant 
from a subset of our ensemble. We subtracted the normal component of the momentum squared (which is typically less than 
10) such that negative values mean refraction and positive ones reflection. For mass ruq = 0.03 on the left we have a number 
of events in the upper left quadrant that would have tunneled with the exact change of the determinant and only a few that 
actually tunneled (in the two lower quadrants). For rnq = 0.05 the picture is similar, even though there are more tunneling 
events. 



For our analysis, we typically discard the first 100 configurations in a stream and separate two consecutive mea- 
surements by 5 trajectories. The separation of the configurations is based on our measurement of the auto-correlation 
time of the plaquette. It is around 5(1) varying little with the quark mass. We observed no significant difference 
in acceptance rate between the quark masses. However, the auto-correlation time of the topological charge differs 
enormously. This is due to the fact that the estimation of the step height of the fermionic action intrinsic to HMC is 
much more subject to fluctuation for smaller masses than for larger ones. As argued in our previous publication, a 
poor estimator results in reflections. We partially address this problem by the use of multiple pseudo-fermion fields, 
but the problem remains. In Fig. ^ we show a scatter plot of the real change in the determinant at the step, from Eq. 
1121 as compared to the stochastic estimate with three pseudo-fermion fields. We subtracted the normal component 
{N,tt)'^ of the momentum so that, according to Eq. negative values allow refractions. However, {N,tt)'^ is almost 
always smaller than 10. The stochastic estimate of the step height has a wide spread but is typically large. Since 
exp(— AS*) will average to the ratio of the fermionic determinants on both sides of the step, this is a consequence of 
the large fiuctuation in the estimator. (A few small values of AS have to be compensated by a large number of large 
ones, for which exp(— AS*) is approximately zero.) 

The low correlation between the estimator and the physical step height Eq. H12|l shows up in the large auto- 
correlation time of the topological charge, whose time history is shown in Fig.|21 Even though part of it is physics — 
lighter quarks make it harder to get from, e.g., = to = ±1 — the height of the step grows with 1/to^ instead 
of the expected determinant ratio, log m. Since the normal component of the momentum is roughly independent 
of the quark mass, it becomes more and more difficult to change topology (also see discussion in Sec. 0). The 
large auto-correlation time for the topology is a phenomenon that is also known with other fermions, e.g. improved 
staggered quarks. To the extent that these formulations know about topology, the step in the fermion action for 
the overlap might be replaced for them by a steep region which approximates the step. The result is the same: if 
the approximation of the determinant is bad, the step is overestimated most of the time and one does not change 
topology. 

Let us finally take a look at the relative cost of the various ingredients of the algorithm. In Table ITTl we list the cost 
per call and the fraction of the total cost of each of the major parts of the program. The conjugate gradient inversion 
of H^^ needed for the computation of the force, starting action and the reflections/refractions takes by far the largest 
fraction. Even though one inverts only on one source (the zero-mode of h(—Ro)) in the reflection/refraction routine, 
these inversions are very expensive since there is no good starting vector. Therefore, they alone take a fifth to a 
quarter of the total cost, depending on the quark mass. The inversions can be cheaper for the lighter quark mass 
because there is less topology; the conditioning number of ifov for v = is lower than for > 0. 

The cost of computing the 8 eigenvectors of the kernel operator is small, about 10%. It is cheap because the 
eigenmodes of the kernel do not change much during the evolution; one thus has good starting vectors. We need 
them to precondition the construction of the overlap operator and to monitor whether an eigenvalue has changed 
sign. Finally, the computation of the fermion force (outside the inversion of H^^) takes about a sixth of the total 
time. This is due to the inversion of the kernel operator, the computation of ipf6h{—Ro)^pi for each of the poles in the 
Zolotarev approximation of the sign function and the inversion of h{—Ro) for the projector term in the sign function 
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FIG. 2: The history of the topological charge u as given by the index theorem for our three quark masses. The data shown 
includes the thermalization runs. 
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sec. per call 


% of total time 


sec. per call 


% of total time 


CG in MD evolution 


406 


48% 


534 


46% 


reflect. / refract. 


800 


18% 


1029 


26% 


fermion force 


148 


17% 


175 


15% 


eigenvalues of h{—Ro) 


110 


12% 


115 


9% 


setup source 


577 


3% 


634 


3% 


gauge update 


27 


3% 


27 


2% 



TABLE II: The cost in wall clock time for various parts of the algorithm measured on a subset of the data. The inversion of 
H^^ takes the largest part. It is responsible for the majority of the cost in the flrst, second and flfth row, together about 75%. 



as discussed at the end of Section Hll 



IV. SETTING THE SCALE 



In order to get an idea at which lattice spacing we are simulating, we measured the heavy quark potential and 
extracted the Sommer parameter vq 25]. (It is defined as the distance at which the force between two heavy quarks 
F{r) with distance r satisfies rQF(ro) = L65.) This is necessary because we expect a substantial shift due to the 
dynamical fermions [2^ with respect to the quenched lattice spacing [2^. Unfortunately, our lattice is relatively 
small. It is thus not possible to estimate the uncertainty in the extraction of the potential from the data alone. To 
illustrate this we show in Fig. O the potential extracted from Wilson loops W^[r, t] which have been constructed from 
HYP smeared links'^sj. The two sets correspond to single exponential fits to the W[r, t] with t g [2, 4] and t g [3, 4], 
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FIG. 3: The heavy quark potential extracted from fits to Wilson loops constructed from HYP links. The two sets result from 
two different fit ranges; t G [2,4] and t £ [3,4]. The curves are fits of the form Eq. 11311 to the potential which include a 
perturbative correction to account for the HYP links. 



(3 = 7.2 


i e [3,4] 


t G [2,4] 


mq 


ro/a 


a[fm] 


ro/a 


a[fm] 


0.03 
0.05 
0.10 


3.27(3) 
3.17(5) 
3.00(4) 


0.153(2) 
0.158(2) 
0.167(2) 


3.04(2) 
2.97(3) 
2.82(3) 


0.165(1) 
0.168(1) 
0.177(2) 



TABLE HI: The Sommer parameter ro and the lattice spacing a in fm from ro = 0.5 fm. We show results from two fit ranges 
used to extract the potential. The fit range for the fit to the potential was r £ [1.4, 6.1] with little variation between different 
choices for this range. 

respectively. The latter set gives a significantly smaller potential. The curves represent fits of the form 

A 

V{r) = - + Br + C + Df{r) (13) 
r 

where /(r) is a perturbatively determined correction to account for the effect of the HYP links on short distances. 
The results are given in Tab. Illll 

To estimate the systematic error from the small volume, we have generated a quenched set of 8^ and 12* lattices 
at f3 = 7.77, Mo = 0.887 with a similar lattice spacing (a = 0.163(1)). Analyzing these lattices, we find that the ro 
extracted from the 8"* lattices with fit range t G [3,4] is the same within error-bars as the one extracted from the 12"* 
and fit ranges starting at 3 or 4 ranging to 5 or 6. In the following, we will therefore work with the lattice spacing 
extracted form the t G [3, 4] fit range. 

The dimensionless quantity TQ^fa has been used in the past to quantify the impact of dynamical quarks on the 
shape of the potential We find ro^/a = 1.10(1) almost independent of the sea quark mass and 1.18(1) for our 
quenched ensemble. On larger lattices and a finer lattice spacing, Ref. [3(1 [ found a quenched value of about 1.16 and 
1.128 for two flavors of dynamical staggered quarks. From Ref. =3]| we get a value of about 1.14 on at a similar lattice 
spacing with two flavors of dynamical clover Wilson fermions. Given the systematic and statistical errors these values 
agree well with our flndings. 

V. TOPOLOGICAL SUSCEPTIBILITY 

The topological susceptibility illustrates the strengths and weaknesses of our simulation. In contrast to all simu- 
lations with non-chiral actions, the measurement of the topological charge in an overlap simulation is trivial. It can 
even be done during the simulation by monitoring zero crossings of the smallest eigenmode of the kernel operator (if 
the topological charge has been determined once at the beginning of the simulation). However, as we have already 
remarked, the autocorrelation time of the topological charge during the simulation is annoyingly long. 

We begin by showing (in Fig. [21 time histories of the topological charge for the different simulations performed at 
our three values of dynamical quark mass. Histograms of the topological charge (including the thermalization runs) 
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FIG. 4: Monte Carlo simulation time between topology changes versus quark mass. 



are shown in Fig.|Sl The topology is recorded at the end of each trajectory (rather than at the end of each time-step). 
No sophisticated analysis is needed to see that the autocorrelation time grows as the quark mass falls, and Fig. 0] 
shows that the mean time between topological changes varies inversely with the square of the quark mass. 

With such long times between tunnelings, we are concerned about thermalization effects in our data. At aruq = 0.1 
{v) is zero within statistical uncertainty (it is -0.05(13) throughout the run) and {v^) seems to have stabilized after 
50 trajectories are discarded, so we cut the data there. At arUq — 0.03 {v) is -0.15(11) when 100 trajectories are 
discarded from each run and 0.04(12) when 150 are discarded; {iP) shows little variation with cuts of more than 50 
initial trajectories per stream, and so we dropped the first 100 trajectories. At aniq = 0.05 the situation is similar; 
we see little variation cutting more than 90-100 configurations and again dropped 100 before averaging. 

We find topological susceptibilities of = 2.17(29) x lO"'' at am, = 0.1, xa" = 1.37(39) x lO"'' at am, = 0.05, 
and xa^ = 1.02(24) x lO""* at am, = 0.03. 

For the quenched data we are able to space configurations used in the analysis far enough apart in simulation time 
as to be essentially uncorrelated. We measured a lattice topological susceptibility of xo-^ = 6.13(76) x lO""*. With the 
quenched rp/a = 3.08, this is x ~ (191 MeV)**, which is quite consistent with typical quenched results, e.g. 32]. 

We attempt to translate this data into dimensionless units in order to facilitate comparisons with other measure- 
ments. We take our measurements of rp/a from the previous section to compute x^o ^-nd do the same with our quark 
masses, using the Z-factor as described in the Appendix to convert them to /i = 2 GeV MS values. We present our 
results in Fig. Diirr|33l| has presented a phenomenological interpolating formula for the mass dependence of the 
topological susceptibility, in terms of the condensate S and quenched topological susceptibility Xq, 

I Nf I 

- = ^ + — . 14 

X Xq 

Taking E from our RMT analysis in the next section ( rg'^S = 0.43) produces the curve shown in the figure. 

Lattice results presented elsewhere typically use the pseudo-scalar mass as the ordinate. As we will describe in the 
next section, we do not have reliable pseudo-scalar masses because our lattices are too small. However, we can use the 
Diirr formula as a benchmark to compare with other simulations. Our results are in rough agreement in magnitude 
with those of another dynamical overlap simulation, at lattice spacing a ~ 0.25 fm, by Fodor, et. al. 9J. The two 
simulations both lie below the Diirr curve. 

Data from simulations with two fiavors of ordinary staggered simulations and was analyzed and published by (among 
others) A. Hasenfratzjs^. At finite lattice spacing all this data lies far above the Diirr curve, and is not too different 
from the quenched result. Other dynamical fermion data with non-chiral fermionsl35| is also high with respect to the 
Diirr formula and to overlap data. It is easy to imagine that simulations with non-chiral actions would overestimate 
the topological susceptibility since their massless Dirac operators do not have exact zero modes. 

We do not think that the small volumes of our simulation have suppressed the susceptibility since our quenched 
simulations are not anomalously low. Little is known about the scaling properties of overlap actions, and our results 
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FIG. 5: Histograms of topological charge from (a) arriq — 0.03, (b) aruq — 0.05, (c) airiq = 0.1 and (d) quenched simulations. 
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FIG. 6: Topological susceptibility versus quark mass, in units of ro. The curved line is the Diirr interpolating formula, Eq. 1141 
The three horizontal lines give the quenched value and its error. 



and those of Ref. |9| have large lattice spacings. We are of course not satisfied with the quahty of our data from the 
point of view of autocorrelations, lattice spacing, extraction of hadron masses, and simulation volume. 



VI. MESON TWO-POINT FUNCTIONS 



Let us now turn to meson two-point functions. The purpose of doing so is two-fold. First, we want to get an idea of 
the pion masses at which we are simulating. This will not work very well since the three volume is small and the time 
extent is far too small for the excited states to decay, but it can provide us with an upper limit of the pseudo-scalar 
mass. The second purpose is to compare the two-point functions to quenched resuhs on similarly sized lattices and 
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FIG. 7: The pseudo-scalar zero momentum two-point function for = 0.03 and niq — 0.05. The lines represent a fit to a 
single cosh in the range t £ [2, 6]. 
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FIG. 8: The squared mass of the pseudo-scalar meson as a function of the bare quark mass. It does not follow the Gell-Mann- 
Oakes-Renner relation. 

look for effects of the dynamical fcrmions in the scalar correlator. We compute zero- momentum correlators 

c^Jit) = :^^(Vi(o,o)r,v(o,o) vi(x,t)r,V'(x,i)) (i5) 

X 

where the fermion fields tp are contracted with the appropriate flavor structure. We compute the quark propagators 
with Gaussian sources of radius 2a on gauge configurations in Coulomb gauge. We use point sinks and apply low-mode 
averaging using the four lowest eigenmodes of the Dirac operator 0, ^| . 

In Fig. Owe show the pseudo-scalar two-point function for our two smaller dynamical quark masses. Unfortunately, 
the corresponding masses, shown in Fig. |S1 do not differ much. This can be attributed to the fact that T = 8 is just 
too small. We thus see a superposition of excited states with the ground state which does not depend on the quark 
mass. 

More impressive is the comparison of the scalar correlator from the dynamical and the quenched ensemble. It is 
known that this correlator turns negative in the quenched theory, which is a sign that quenched QCD is not unitary. In 
Fig. (bottom) we show that this is the case even on a small lattice. The full two-point function is shown (connected 
by a line) as well as the contributions of the various topological sectors. Whereas the function is positive in i/ = 0, the 
zero-modes in the sectors of non-trivial topology turn it negative. In the full theory, however, the fermion determinant 
suppresses the sectors with \i>\ > sufficiently for the whole two-point function to remain positive (Fig.|51upper row). 

VII. THE CONDENSATE FROM RANDOM MATRIX THEORY 



It was proposed more than a decade ago that the distribution of the low-lying eigenvalues of the QCD Dirac operator 
in a finite volume can be predicted by random matrix theory (RMT) |38l l39l |40j . Since then this hypothesis has 
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FIG. 9: The scalar zero momentum two-point function for the dynamical ruq = 0.03 and — 0.05 ensembles and for the 
quenched ensemble with a valence quark mass of rUq — 0.03. Contrary to the quenched theory, the scalar two-point function 
stays positive with dynamical quarks. We show the full two-point function (connected by a line) and the contributions from 
the each topological sector separately. The topologically non-trivial configurations turn the correlator negative in the quenched 
theory, where they are more abundant. 



received impressive support from lattice calculations, mainly quenched simulations |4lL l42l l43l l44l l45l | . but also some 
dynamical ones using staggered quarks 0, 0| . 

Typically, the predictions are made in the so-called epsilon regime, for which 1/A <C i <C I/tItt with A a typical 
hadronic scale. However, it has been found that they describe the data in a wider range. In a recent large scale study, 
e.g., using the overlap operator on quenched configurations [i^, it could be shown that from lattices with a length 
larger than 1.5 fm, the RMT predictions match the result of the simulation. Our dynamical lattices have a spatial 
extent of about 1.3 fm. As we will see, random matrix theory describes our low- lying Dirac spectra quite well. 

Our analysis is based on the distribution of the k-th eigenmode from RMT as presented in Ref. |3 and successfully 
compared to simulation results in Ref. [4^ . The prediction is for the distribution of the dimensionless quantity 
( — AfcSy in each topological sector with Xk is the fc-th eigenvalue of the Dirac operator, E the chiral condensate and 
V the volume of the box. These distribution are universal and do not depend on additional parameters other than 
the number of flavors, the topological charge and the dimensionless quantity mnSV. (Note that topology effects the 
distributions, in contrast to the behavior of staggered fermions seen by Ref. |47l.) By comparing the distribution of 
the eigenmodes with the RMT prediction one can thus measure the chiral condensate S. The main advantage of this 
method is that it gives the zero quark mass, infinite volume condensate directly. The validity of the approach can be 
verified comparing the shape of the distribution for the various modes and topological sectors. The main uncertainty 
comes from a too small volume which causes deviations in the shape, particularly for the higher modes. As we will 
see below, the direct measurement of S from {ip'ip)(m) ^ which requires extrapolation into the chiral limit as well as 
correction for the finite volume, is unreliable for our data. 

In Fig. ^] we show the distribution of the two lowest eigenmodes of the overlap operator (scaled by Ty) measured on 
the V = Q and v = ±\ parts of the aniq = 0.03 and aruq = 0.05 ensembles. To correct for the fact that our eigenvalues 
lie on a circle of radius Rq instead of a straight line, we use the stereographic projection A = | A| \Jl — |Ap/2i?o. We fit 
the RMT prediction from Ref. to these distributions. The chiral condensate S is the only free parameter in this fit; 
we get Ey/a = 54(2) and EF/a = 53(2), respectively. This corresponds to rgE = 0.46(2) and 0.41(2). We studied the 
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dependence of this result on the number of trajectories used in the thermaHzation and the separation of consecutive 
configurations. We found no systematic variation beyond omitting the first 100 trajectories and separating them by 5. 
The prediction agrees overall well with the measured distribution given the low statistics. However, the distribution 
of the lowest mode in the v — sector seems to have a tail at larger X'EV that does not match the prediction. This 
could be an effect of the small volume. We also show the prediction for the distribution of the third mode from our 
fitted values of SV^ in the third column of Fig. ^| The RMT curve and the data, again, agree quite well. However 
for the = 1 sector, the curve seems to be on the right of the data. This is probably a sign of the break-down of 
RMT for eigenvalues larger than the Thouless energy 43, 49]. 

The RMT predictions are made with the assumption that the volume is infinite. We (obviously) are not in that 
situation. In finite volume, in the epsilon-regime of chiral perturbation theory, finite volume modifies the formula for 
the condensate by multiplication by a shape factor, E pE, where 

and c{li/l) depends on the geometry jS^. We do not know p since we have not measured /^, but combining our lattice 
spacing and lattice size with = 93 MeV gives p ^ 1.4. This is too large a correction to be trustworthy; again, we 
need a lattice with a larger physical size. Our two rgE values must be multiplied by Ss and p to give a continuum 
value: with tq = 0.5 fm, pE(MS) = 0.032(1) or 0.029(1) GeV^ or {pj:(MS)y/^ ^ 316 or 304 MeV. (Errors are from 
To/a and the fit uncertainty in S.) 

We also attempted to measure the mass dependent quark condensate {'iptpim)} using 12 random sources per lattice 
on subsets consisting of about 50 lattices per mass value, spaced by 5 trajectories. Our data is displayed in Fig. ^2 We 
show the mass-dependent condensate summed over all topological sectors as well as its value in sectors of topological 
charge ly — and — 1 (with the zero mode contribution removed for the latter case). As expected, there is no 
divergence in {ipip{m)) at zero quark mass (as would occur in a quenched simulation). 

The parameter S must be extracted from a fit of {tpi'ijn)) to some functional form. For example, were we in a 
large volume, we would fit to the usual constant plus chiral logarithm plus Crriq formula. 

We believe that we are, in fact, in a small volume, and that the appropriate fitting function is that for the so-called 
e— regime, given by Gasser and Leutwyler^ for the average-topology case, or by Leutwyler and SmilgalSJ, for the 
case of fixed topology. In either case condensate is given by 

, p^dZ 

with Z the appropriate partition function written by the above authors or given equivalently by random matrix theory, 
and fi = mqT,V . We append a linear Cruq term to Eq. 1171 and attempt to fit the parameters S and C to the data. 
This was not successful. 

First, we do not know which of our mass values lie in the region of validity of our fitting function. Fig. 1121 shows 
a set of fits to the v — condensate. The fitting function is a sum of the fixed-topology expression plus a linear 
mass-dependent term. Curves (a) and (c) show the result of fits to the three or lowest two mass data points. Curves 
(b) and (d) are the part of the fit coming from the S(y, /i) expression of Eq. El The three-mass fit gives rgS = 0.53(4) 
while the two-mass fit gives 0.67(8). These numbers are obviously unstable. We are currently generating data at 
smaller quark masses which should alleviate this problem by filling in the curve. 

But the numbers we record are quite different from the RMT ones. We think (see Ref. [52| for a useful discussion) 
that this is an artifact of the small simulation volume. The problem of principle that we face is that the formulas to 
which we fit (V'?/'(m)) assume that the condensate is given by an integral over a known spectral density p^^^Q, p) for 
all values of C, 



2W dC^^:^. (18) 



We are doing simulations in small coarse lattices, and so this assumption is unwarranted. We expect that cutoff effects 
will alter the high eigenvalue part of the spectrum. Only a larger lattice will cure this problem. 

We can check this assumption by fitting (■^■^(to)) in the \v\ = 1 sector. We get E = 0.70(3) from a fit to the 
am — 0.03, 0.05 and 0.1 data sets, and 0.88(10) from the am — 0.03 and 0.05 sets. These are bigger discrepancies 
from the RMT results than the v = Q fits. Fig. IIUI shows that the third eigenmode of the v = Q sector matches the 
RMT prediction much better than the third \h'\ = 1 mode does. 
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FIG. 10: The distribution of the lowest two eigcnmodcs of the Dirac operator for our ensemble for the sector of trivial topology 
and V = ±1. The lines are the result of fits of the random matrix theory prediction to the data for the two lowest modes. They 
correspond to T,V/a = 54(2) and T,V/a = 53(2) for the amq = 0.03 and aruq = 0.05 ensemble, respectively. The lines for the 
third mode are predictions. 



VIII. CONCLUSION 



We have presented results from simulations of two flavors of dynamical overlap fermions at sea quark masses of 
about 35, 55 and 100 MeV. For us, this is a second step in gaining experience with these simulations. We therefore 
focused on dynamical fcrmion effects in the results. We were able to show that the topological susceptibility is greatly 
reduced as compared to quenched simulations. This measurement also is a great strength of the use of overlap fermions 
for the sea quarks. The topological charge as defined by the index theorem has a direct impact on the update of 



14 



i 
<> 



o 

I 

5 



0.0 0.1 0.2 0.3 0.4 

FIG. 11: The mass- dependent condensate vs quark mass, scaled by appropriate powers of ro, from direct measurement. 
Octagons show the full result (summed over all topological sectors). Squares and diamonds show tp^p{m)) in sectors of topological 
charge |;^| = and 1 respectively, the latter with the zero-mode contribution removed. 
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FIG. 12: Examples of fits to the v = mass-dependent chiral condensate. Curves (a) and (c) use the three mass points and 
the two lower ones, respectively. Curves (b) and (d) show the part of the fit coming from the S( V, jj,) expression of Eq. 1171 

the gauge fields during the simulation. However, we are still worried about the long auto-correlation time of the 
topological charge. 

We also extracted the quark condensate from a comparison of the distribution of the lowest eigenvalues with random 
matrix theory. This method is much simpler than a direct fit to {ipip(rn)) and has the advantage of giving the infinite 
volume, zero mass condensate directly without need of extrapolation. The eigenvalue distributions fit consistently the 
distributions in the various topological sectors for the first three modes. However, due to our very small volume, there 
is a finite size correction to E which is not well under control. This correction is probably 0(40%) for our simulation. 
But it is expected to scale which 1/L^ and one thus only needs a moderately larger volume to make it small. 

We were also able to demonstrate the effect of the dynamical quarks on the scalar meson two-point function. We 
showed that contrary to the quenched theory the scalar two-point function is not turned negative by the zero-modes 
in the sectors of non-trivial topology. 

All this is possible through an efficient implementation of the overlap operator using fat links, multiple pseudo- 
fermion fields to decrease the auto-correlation time of the topological charge and improvements in the algorithm 
namely the computation of the height of the step in the fermionic action. 
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process 


Zi 


q*a 


Zv,A 

Zp,s 


0.54 
2.96 


2.24 
2.85 



TABLE IV: Table of Z-factors and (j*'s for our action, defined so Zi = 1 + Zig^ {q*)Cf /{16tt^) = 1 + Zia(g*)/(37r). 



APPENDIX A: Z-FACTORS FROM PERTURBATION THEORY 

A simple application of the techniques described in Ref. [s^l gives us one loop predictions for the vector, axial 
vector, pseudo-scalar, and scalar currents for our action. They are shown in Table HVI The value of the momentum 
scale at which the strong coupling constant is evaluated (from the Lepage-Mackenzie convention 54J ) is also shown. At 
one-loop order there is no difference between the gauge propagator from a tadpole improved action and the tree-level 
one, so we use the tree-level propagator in our computation. In this order of perturbation theory, each step of stout 
smearing is equivalent to a step of unitarized APE-smearing|55j with p = a/6 in the terminology of Ref. [s^ . 

In practice, we define the coupling through the so-called "ay" scheme. We only know the one-loop expression 
relating the plaquette to the coupling; it is 

\n^TrUp^-^av{q*)W (Al) 

with W = 0.366 and q*a — 3.32 for the tree-level Liischer Weisz action. In our calculation of the condensate, we 
take the lattice spacing from the Sommer parameter. The coupling from the plaquette is matched to its MS value 
and run to the needed value of q*, where we perform the match. Then the MS result is run to ^ = 2 GeV using 
the usual two-loop formula. In our simulations ay(3.32/a) — 0.192, 0.193, 0.193, and Zg = 1.19, 1.22 and 1.23 for 
the amq = 0.03, 0.05 and 0.10 data sets. Essentially all the difference from unity comes from the MS running from 
12 = 1/a to 2 GeV. Using — ^/Zs gives the MS quark masses quoted in the body of the paper. 



ACKNOWLEDGMENTS 



This work was supported by the US Department of Energy. It is a pleasure to thank A. Hasenfratz for conversations 
and to P. Damgaard for extensive correspondence about random matrix theory. Simulations were performed on the 
University of Colorado Beowulf cluster. We are grateful to D. Johnson for its construction and maintenance, and to 
D. Holmgren for advice about its design. 



P. H. Ginsparg and K. G. Wilson, Phys. Rev. D 25, 2649 (1982). 

H. Neuberger, Phys. Lett. B 417, 141 (1998) arXiv:hep-Iat/9707022l 

H. Neuberger, Phys. Rev. Lett. 81, 4060 (1998) arXiv:hep-Iat/9806025'. 

P. Hasenfratz and F. Niedermayer, Nucl. Phys. B 414, 785 (1994) a rXiv:hep-lat/9308004^ . 

D. B. Kaplan, Phys. Lett. B 288, 342 (1992) arXiv:hep-lat/920601 3 | . " 

V. Furman and Y. Shamir, Nucl. Phys. B 439, 54 (1995) arXiv:hep-lat/9405004l. 

A. Bode, U. M. Heller, R. G. Edwards and R. Narayanan, arXiv:hep-lat/9912043 

Z. Fodor, S. D. Katz and K. K. Szabo, JHEP 0408, 003 (2004) arXiv:hep-lat/03 lToTo] . 

Z. Fodor, S. D. Katz and K. K. Szabo, arXiv:hep-lat/0409070 

T. DeGrand and S. Schaefer, Phys. Rev. D 71, 034507 (2005) arXiv:hep-lat/0412005 . 

N. Gundy, S. Krieg, G. Arnold, A. Frommer, T. Lippert and K. Schilling, arXiv:hep-lat /0502007| 

C. Aubin et al. [MILC Collaboration], arXiv:hep-lat/0409051 

C. Bernard et al., Phys. Rev. D 68, 114501 (2003) arXiv:hep-lat/0308019'. 

C. Morningstar and M. J. Peardon, Phys. Rev. D 69, 054501 (2004) arXiv:hep-Iat/0311018]. 

S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth, Phys. Lett. B 195, 216 (1987). 

M. Hasenbusch, Phys. Lett. B 519, 177 (2001) arXiv:hep-lat/0107019 . 

M. Hasenbusch and K. Jansen, Nucl. Phys. B 659, 299 (2003) arXiv:hep-lat/0211042 . 

P. Hasenfratz, V. Laliena and F. Niedermayer, Phys. Lett. B 427, 125 (1998) arXiv:hep-lat/9801021|. 

G.H. Golub and C.F. van Loan, Matrix Computations, 3rd ed., Baltimore: The Johns Hopkins University Press, 1996 

S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken and R. L. Sugar, Phys. Rev. D 35, 2531 (1987). 

R. Narayanan and H. Neuberger, Phys. Rev. D 62, 074504 (2000) arXiv:hep-lat/0005004 . 



16 



M. Liischer and P. Weisz, Commun. Math. Phys. 97, 59 (1985) [Erratum-ibid. 98, 433 (1985)]. 

M. G. Alford, W. Dimm, G. P. Lepage, G. Hockney and P. B. Mackenzie, Phys. Lett. B 361, 87 (1995) 
arXiv:hep-lat/9507010 . 

J. C. Sexton and D. H. Weingarten, NucL Phys. B 380, 665 (1992). 

R. Sommer, NucL Phys. B 411, 839 (1994) arXiv:hep-lat/93 10022 . 

A. Hasenfratz and T. A. DeGrand, Phys. Rev. D 49, 466 (1994) arXiv:hep-lat /9304001" . 
C. Gattringer, R. Hoffmann and S. Schaefer, Phys. Rev. D 65, 094503 (2002) arXiv:hep-lat /Ol 12024] . 

A. Hasenfratz and F. Knechth, Phys. Rev. D 64, 034504 (2001) arXiv:hep-lat/0103029 . 

A. Hasenfratz, R. Hoffmann and F. Knechth, NucL Phys. Proc. SuppL 106, 418 (2002) .arXiv:hep-lat /Ol 10168) . 
C. W. Bernard et ai, Phys. Rev. D 64, 054506 (2001) arXiv:hep-lat/0104002 . 

A. Ah Khan et al. [CP-PACS CoUaboration] , Phys. Rev. D 65, 054505 (2002) [Erratum-ibid. D 67, 059901 (2003)] 
arXiv:hep-lat/0105015 . 

C. Gattringer, R. Hoffmann and S. Schaefer, Phys. Lett. B 535, 358 (2002) |arXiv:hep-lat/0203013| . 
S. Diirr, NucL Phys. B 611, 281 (2001) arXiv:hep-lat/0103011 . 
A. Hasenfratz, Phys. Rev. D 64, 074503 (2001) arXiv:hep-lat/0104015 . 

See, for example, C. R. AUton et al. [UKQCD Cohaboration] , Phys. Rev. D 70, 014501 (20 04) |arXiv:hep-lat/0403007| . 
T. DeGrand and S. Schaefer, Comput. Phys. Commun. 159, 185 (2004) arXiv:hep-lat/04010m 
T. DeGrand and S. Schaefer, NucL Phys. Proc. SuppL 140, 296 (2005) arXiv:hep-lat/0409056 . 
E. V. Shuryak and J. J. M. Verbaarschot, NucL Phys. A 560, 306 (1993) arXiv:hep-th/9212088|. 
J. J. M. Verbaarschot and I. Zahed, Phys. Rev. Lett. 70, 3852 (1993) arXiv:hep-th/9303012j: 
J. J. M. Verbaarschot, Phys. Rev. Lett. 72, 2531 (1994) arXiv:hep-th/9401059 . 

M. E. Berbenni-Bitsch, S. Meyer, A. Schafer, J. J. M. Verbaarschot and T. Wettig, Phys. Rev. Lett. 80, 1146 (1998) 
arXiv:hep-lat/9704018 . 

P. H. Damgaard, U. M. Heher and A. Krasnitz, Phys. Lett. B 445, 366 (1999) arXiv:hep-lat/9810060|. 
M. Gockeler, H. Hehl, P. E. L. Rakow, A. Schafer and T. Wettig, Phys. Rev. D 59, 094503 (1999) arXiv:hep-lat/9811018l. 
R. G. Edwards, U. M. Heher, J. E. Kiskis and R. Narayanan, Phys. Rev. Lett. 82, 4188 (1999) arXiv:hep-th/9902117^ . 

L. Giusti, M. Liischer, P. Weisz and H. Wittig, JHEP 0311, 023 (2003) arXiv:hep-lat/0309189 

M. E. Berbenni-Bitsch, S. Meyer and T. Wettig, Phys. Rev. D 58, 071502 (1998) arXiv:hep-lat /980 4030| . 

P. H. Damgaard, U. M. Heller, R. Niclasen and K. Rummukainen, Phys. Lett. B 495, 263 (2000) ar Xiv:hep-lat/0007041| . 
P. H. Damgaard and S. M. Nishigaki, Phys. Rev. D 63, 045012 (2001) arXiv:hep-th/0006111 . 
J. C. Osborn and J. J. M. Verbaarschot, Phys. Rev. Lett. 81, 268 (1998) arXiv:hep-ph/9807490l 

J. Gasser and H. Leutwyler, Phys. Lett. B 184, 83 (1987); NucL Phys. B 307, 763 (1988). See also P. Hasenfratz and 
H. Leutwyler, NucL Phys. B 343, 241 (1990). 
H. Leutwyler and A. Smilga, Phys. Rev. D 46, 5607 (1992). 

For a longer discussion of these issues, see P. H. Damgaard, NucL Phys. Proc. Suppl. 106, 29 (2002) |arXiv:hep-lat/0110192| . 
T. DeGrand, Phys. Rev. D 67, 014507 (2003) arXiv:hep-lat/0210028 . 
G. P. Lepage and P. B. Mackenzie, Phys. Rev. D 48, 2250 (1993) arXiv;hep-lat/9209022 . 

M. Albanese et al. [APE Collaboration], Phys. Lett. B192, 163 (1987); M. Falcioni, M. L. Paciello, G. Parisi and B. Tagli- 
enti, NucL Phys. B251 (1985) 624. 
[56] C. W. Bernard and T. DeGrand, NucL Phys. Proc. Suppl. 83, 845 (2000) .arXiv:hep-lat/9909083| . 



